Alberta Thy 02-99 
Even-Odd and Super-Even Effects in the Attractive Hubbard Model 



K. Tanaka and F. Marsiglio 
Department of Physics, University of Alberta, Edmonton, Alberta, Canada T6G 2J1 

(February 1, 2008) 






O 
o 

Oh 



X3 
O 



> 

en 

(N 

O 
0^ 
On 
-I— > 



I 

O 

o 



X 



The canonical BCS wave function is tested for the attrac- 
tive Hubbard model. Results are presented for one dimen- 
sion, and are compared with the exact solutions by the Bethe 
ansatz and the results from the conventional grand canonical 
BCS approximation, for various chain lengths, electron den- 
sities, and coupling strengths. While the exact ground state 
energies are reproduced very well both by the canonical and 
grand canonical BCS approximations, the canonical method 
significantly improves the energy gaps for small systems and 
weak coupling. The "parity" effect due to the number of elec- 
trons being even or odd naturally emerges in our canonical re- 
sults. Furthermore, we find a "super-even" effect: the energy 
gap oscillates as a function of even electron number, depend- 
ing on whether the number of electrons is Am or 4m -I- 2 (m 
integer). Such oscillations as a function of electron number 
should be observable with tunneling measurements in ultra- 
small metallic grains. 

PACS number(s): 74.20.Fg, 71.24.-|-q, 71.10.Fd, 71.10.Li 



I. INTRODUCTION 

The possibility of fabricating ||l[ tunnel junctions con- 
taining nanoscale particles has opened up new areas of 
exploration. Experiments can now probe the energy 
spectrum in these particles and study changes as a func- 
tion of temperature and magnetic field |^,|| . In this way 
one can observe changes in the spectrum that are ex- 
pected to occur in the bulk due to phase transitions, and 
determine to what extent these concepts are relevant for 
small systems H] . The unprecedented control over exper- 
imental conditions — changes due to a single tunneling 
electron can be observed — allows one to ask and address 
questions which, until now, have been academic only. 

In this study we follow the experiments of Tinkham et 
al. ^^ and address the issues concerning superconduc- 
tivity in small systems. Specifically they studied small 
superconducting Al particles (diameter of order 10 nm) 
and probed, through tunneling experiments, the excita- 
tion spectrum as a function of temperature and magnetic 
field. In attempting to treat small electronic systems, 
there are many concerns related to possible surface and 
impurity effects, both of which could lead to localiza- 
tion, for example. For the present we ignore these poten- 
tial complications, and instead focus on a question which 
arises in the application of the Bardeen-Cooper-Schrieffer 
(BCS) jsl theory of superconductivity: to what extent is 
the grand canonical ensemble useful (which in the ther- 



modynamic limit, is equivalent to the canonical one) for 
systems with a small number of electrons? 

We propose to tackle this question in a systematic way, 
using the attractive Hubbard model. The choice of this 
model is motivated by several factors. It has long served 
as a paradigm for s-wave superconductivity, and serves 
as the 'minimal' model that best describes superconduc- 
tivity. All the energy scales are very well defined in the 
problem, so that no high frequency (and ill-defined) cut- 
offs are required g] . Exact solutions are available in one 
dimension via Bethe Ansatz techniques (7| , and the grand 
canonical BCS solutions have recently been evaluated for 
large system sizes [g| . In fact in that work some prelimi- 
nary canonical solutions were also examined, but only for 
very small system sizes. In this work we reformulate the 
canonical solution, in such a way that larger systems up 
to three dimensions can be tackled. We note that after 
this work was completed, a paper appeared M in which 
the canonical BCS equations were formulated and solved 
for a model system with uniformly spaced level spacings. 
We will comment on the similarities and differences with 
our own work in the course of our discussion below. 

The outline of the paper is as follows. In the next sec- 
tion we formulate the problem. We mention the work of 
Falicov and Proetto pM, who applied a canonical BCS 
method to the problem of electrons on a tetrahedron (a 
'small' fee lattice). We have generalized this method for 
larger systems, and solved the ensuing equations numer- 
ically. However, memory requirements become very de- 
manding, and therefore, beyond a certain lattice size, we 
had to abandon this approach. We include it in the ap- 
pendix, nonetheless, because the variational wave func- 
tion in this method is more general than the BCS canon- 
ical wave function, and the variational equation is linear. 
In practice, however, as will be shown, this wave function 
only marginally improves the ground state energy. 

We then formulate the canonical variational problem 
strictly in terms of the 5(k)'s familiar from the grand 
canonical solution [[ll| . The resulting equations are te- 
dious, and can be cast in a more enlightening form by fol- 
lowing the pioneering work of Dietrich, Mang and Pradal 
|l2[ , who solved the canonical BCS problem for nuclei. 
This "method of residues" is based on representing ex- 
pectation values in the BCS canonical ensemble through 
contour integrals of the BCS grand canonical ensemble 
expectation values: the latter seems to have started with 
a paper by Bayman [n3. This is also the methodology 
adopted by Braun and von Delft [BJ, who numerically 
evaluated residue integrals by fast Fourier transform. As 
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will be explained, however, we chose to evaluate the re- 
quired integrals analytically (which requires numerical 
summations) . A concise summary of the residue method 
and further references are available in Ref. [Oj. 

Finally, we provide a review of the grand canoni- 
cal formulation and the Bethe Ansatz solutions in one 
dimension [[7|p^J^, for purposes of comparison. We 
have also performed some numerical diagonalizations, al- 
though these are limited to small systems and (in one 
dimension) serve only to verify the Bethe Ansatz solu- 
tions. In both the grand canonical and canonical BCS 
formulations we have performed the variation strictly in 
terms of the single set of parameters, ^(k). This departs 
from more conventional treatments where two sets of pa- 
rameters, M(k), and w(k), are used, with an auxiliary 
relation between them. Our formulation is conceptually 
clearer, but slightly more cumbersome, so some details 
are provided. 

The following section is devoted to numerical results. 
All the computational results presented in this paper will 
be in one dimension. This is because the exact results 
for 'large' systems are only available in 1-D, and there- 
fore it is presently only in this case that one can study 
the full crossover from the bulk limit to small systems. 
Because the attractive Hubbard model is a very local 
model, many of the features of the variational solution 
will remain in higher dimension; in fact we expect the 
agreement of the canonical BCS results with the exact 
ones to improve, but we have no way at present to check 
this. In addition, many of the features of the solutions 
in higher dimensions (such as long range order) will not 
be present in one dimension. A more systematic study of 
the canonical BCS solution in two and three dimensions 
will be presented elsewhere (unfortunately without exact 
checks, except for very small systems). 

Aside from evaluations of the canonical BCS formu- 
lation we will demonstrate the so-called 'parity effect' 
|l6| , |l7| , where qualitative differences in the tunneling gap 
can be observed, depending on whether an even or odd 
number of electrons occupies the small Al particle. While 
this has been understood within a parity-conserved grand 
canonical BCS formulation jlj], it becomes much clearer 
in a canonical formulation. In addition we have found an 
interesting Am vs. Am -f 2 effect, where m is an integer. 
In our case we have observed oscillations as a function 
of electron number. Therefore we have variations in the 
gap between odd, even and super-even (i.e. multiples of 
4) electrons. Hints of such behaviour, noted as a function 
of lattice size, were first discussed by Fye et al. p9| . 

The concluding section summarizes our results. In the 
Appendix we outline the formulation of the linearized 
canonical variation. 



II. FORMULATION 



A. Model 



The attractive Hubbard Hamiltonian is given by 
H ^~J2*s (aI+5,o-«»'T + h.c.) -\U\Y^ n.^n^i (la) 



\U 



5Z^kaL«k<T ~ ^Yj 4T«-k+qi«-k'+qiak'T , (lb) 



kk'q 



where a]^ (oio-) creates (annihilates) an electron with spin 
a at site i and riia- is the number operator for an electron 
with spin a at site i {i is the index for the primitive vec- 
tor Ri). The ts is the hopping rate of electrons from one 
site to a neighbouring site R^ away (often nearest neigh- 
bours only are included), and \U\is the coupling strength 
between electrons on the same site; here the fact that it 
is attractive is explicitly included from the beginning. In 
Eq. dig), we have Fourier transformed the Hamiltonian 
with periodic boundary conditions for N sites in each 
dimension. The aj^^ and ak<T are the creation and anni- 
hilation operators in the reciprocal space, and the kinetic 
energy is given by 



Ek = —2 2_, ^i ^'^^ ^ ' -f^* 



(2) 



where R^ is the coordinate vector connecting sites i to 
i + 5. 

We perforin variational calculations using the compo- 
nent of the BCS wave function that has a given number 
of pairs v and thus conserves the number of electrons 
A^e = 2t/. The BCS wave function is a superposition of 
pair states with all the possible numbers of pairs {v}: 



|BCS) 



GC 



n 



Mk 



^kakT«-ki) |0> 



El*2.), (3) 



v=0 



where |0) denotes the vacuum state, and I^I^o) = |0). The 
j^-pair component 1^2!^) can be obtained by rearranging 
|BCS)gc into a power series of the pair creation operator 



jltft n and can be written as 



1*2.) = E E • • • E n (.9(k.)«LT«Ui) io) (4) 

ki < k2 < <ki, i=l 

=^n(E5(k.)<T«-k.i)io) 



i=i ki 

(ki ^ kj , for all i ^ j) 



(4') 



1/^ 



Here we have defined g(ki) = (Ilk'^k) ^ki/uk;- The 
condition ki < k2 < ■ • ■ < k^, for the {ki} sums in 
Eq. (0) means that all the k^'s should be different and 
ordered so that all the different combinations of {ki} are 
counted once each: the latter condition is removed for 
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the sums in Eq. (Q) and the division by v\ compensates 
for the multiple counting. In a finite system with N sites 
in each dimension, the wave number can take N^ values, 
where D is the dimension of the system, and hence the 



v sums over k's in Eq. 
terms. 



result in f^oC, 



N"i. 



B. Grand Canonical Variation 

We begin by reviewing the grand canonical formulation 
of BCS theory. We have opted, contrary to convention, to 
formulate the minimization problem in terms of only the 
variational parameters, gk = 5(k) (and not the Uk's and 
Wk's). This was originally motivated by the loss of clear 
meaning for the Uk's and iik's in the canonical ensemble, 
though in fact one can proceed equally well with these 
121,01. One therefore simply writes: 



|BCS) 



GC 



n(i 



.9kalt"-ki 



|0), 



(5) 



where the coefficient c properly normalizes the wave func- 
tion. One then calculates the various expectation values 
to determine the ground state energy: 



E^ 






(6) 



where, in this case, the wave function is the BCS grand 
canonical one given in Eq. (^. For example, the expec- 
tation value of the kinetic energy operator, T is 



(vI'|r|vI/) = 2|cp[](l+4)^(6k-M) 



9l 



^ + 9t 



(7) 



Note that normally the product and the denominator are 
absent because these factors are normally defined to be 
unity due to normalization. Also note that the gk's are 
generally complex, but can be chosen to be real; this is 
presumed in Eq. (M), and will be implicit in what follows. 
Careful evaluation of the potential energy terms yields 
the following result: 



E, 



GC= 2 ^ (ek - Ai) Y 



9l \U\ 



El 



5k 



9l 



\U\ y^ .9k .gic' \U\ Y^ Sk ffk 



k/k' 



5^1 + ffk' N 



2 2 

„l + 9ll + 9l ^' 



k^k' 



with fi the chemical potential; it plays the mathematical 
role of the Lagrange multiplier for the condition that the 
average electron number Gc(BCS|7Ve|BCS)GC = Ng. We 
have written the total energy term by term in the fol- 
lowing sequence: the first term is the kinetic energy, the 
following two come from Cooper pair scattering (q = 0), 
and the fourth term is the Hartree term (q 7^ 0) . When 



displayed this way it is apparent that the last term ex- 
cludes the spurious Hartree term whereby an electron in- 
teracts with itself (a '1/N effect'). It is also apparent that 
the reduced JiCS Hamiltonian (q = scattering only) will 
exclude this last term, and in doing so omits the Hartree 
term (since it is often deemed to merely represent a shift 
in energies). At the same time Eq. (H) appears to have 
terms of order 1/N smaller than the dominant terms. In 
fact it docs not, and we can rewrite this equation in the 
following form: 



Egc^ 2 ^ (ek - /i) y 



9l 



9l 



M V 

N 2^ 



k,k' 



5k 



ffk' 



9l 



9l 



l + 5kl + 5^' l + 5^1 + .9k' 



(9) 



where now the summations are unrestricted. The en- 
tire expression is now clearly extensive, as it should be, 
though in this form pairing terms of order 1/A'^ have been 
mixed in with the Hartree term. 

The next step is to carry out the variation with respect 
to the gk's; this can be done straightforwardly to yield: 



2(^k-A)..k = ^ETfk[i-^^] 



9k' 



where 



fi = fj,+ 



\U\ Y- 9l 
N ^l+9l 

We also define the pair potential, 

\U\ ^_gk 

k 



A 



BCS 



1^1 y^ 



9l 



(10) 



(11) 



(12) 



so that the BCS equation can be written 

2 (ek - m) .9k = Abcs [l - 5k] ■ (13) 

The solution is 

Ek - (ck - f-t) 



9k = 



iBCS 



(14) 



where E]^ = ^{e^ — ji)'^ + A|(-,g is the quasi-particle en- 
ergy. This solution is only implicit, since Abcs depends 
on the gk's through Eq. dl2), and must be determined 
by numerical iteration. The number equation is also re- 
quired to determine the chemical potential as a function 
of coupling strength. It is |ll| 



^ ^ (^e) ^ ^ I y (£k - fi) 



N 



N ^ Eu 

k 



The gap is then defined to be 

Ao = min (i?k) 



(15) 



(16) 
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If we define the minimum band energy, e^in ^ —D/2, 
with D the bandwidth, then it follows that Ao = 
•\/(emin ^ P-)^ + ^Bcs '^^en jl < eniin- This definition 
implies that the gap depends on a certain momentum, 
determined by the wave vector k at which the minimum 
energy occurs. The total energy also follows readily from 
these expressions. It is 



N 



1 



^E^k 



TV 



Sk 



\U\l- 



IS? 



\u\ 



where we have used the fact that 



9k 



- N EkT+^' 



(17) 



The 



BCS grand canonical results in the remainder of the pa- 
per have been obtained through the solution of these 
equations on finite lattices. 



C. Canonical Variation 

1. BCS 

Rather than adopt the wave function given by Eq. 
(a), which spans all (even) number sectors, one can 
work directly with the wave function |^2i/} given by 
Eq. (Q) with fixed electron number, A'e = ^v. A 'sim- 
plification' occurs if we try to linearize the problem, 
by defining a variational parameter C(ki,k2, ■ ■ • ,k,y) = 
(7(ki) (7(k2) • • ■ .g(kiy) for this wave function, following 
Refs. [|0| and |^ for the two-pair case. Thus one can 
rewrite Eq. (Q) as 

If 
1*2.) = E E • ■ • E ^(ki, k2, • • • , k,) n 4.T«Lka io) 

ki < k2 < < k„ 8=1 

(18) 

As discussed in Ref. Q, the C's defined as above are 
subject to constraints; for example, for i^ = 2 and for 
a given set of four k values (qi, q2, Qs, q4) which satisfy 

qi < q2 < qs < q4, 

C(qi,q2)C(q3,q4) = C(qi, qa) C(q2, q4) 

= C(qi,q4)C(q2,q3). (19) 

If these constraints are ignored and the C's in Eq. (|l8) 
are treated as p^dCu independent parameters, the mini- 
mization of the energy 



-E2 



(^2.|g|^2.} 
(*2.|^2.) 



(20) 



with respect to {C} reduces to a linear problem pOl. 
We initially adopted this method of canonical variation, 
which actually allows more variational freedom than us- 
ing the 5k's. However, as one can readily appreciate. 



the number of variational parameters grows very quickly 
with increasing lattice size, and therefore in practice, 
this method has limited usage. Where it was practical, 
however, we carried out these calculations and compared 
them with the canonical BCS results to be described be- 
low. The gain in accuracy for the ground state energy 
was fairly small. Our formulation and results in this case 
are summarized in the Appendix. 

The cost of working with the gk's is that the problem 
remains a very nonlinear one. Nonetheless the gain in 
computational ease more than compensates for this, in 
that only N variational parameters are required, where 
N is the number of lattice sites. As we shall see, this 
allows us to easily study the limit from small systems 
to those that are well described by the grand canonical 
ensemble. The formulation of this problem using the ex- 
pression for |^2i/) given in Eq. (|j) is straightforward but 
tedious. We therefore follow the "method of residues" 
|12| (see also ||lj] and more recently 1^) originally devel- 
oped for the nuclear pairing problem. An advantage of 
this method is that matrix elements are more easily eval- 
uated (they remain as simple as the ones in the grand 
canonical formulation). Moreover the energy and the re- 
sulting variational equation can be written in a compact 
way which parallels the grand canonical BCS equations. 
The starting point is to write Eq. (0) in the general form 
as a particle-number projection of the grand canonical 
BCS state [O]. In more general terms, this is the pro- 
jection that restores the symmetry of the Hamiltonian 
(i.e., conserved particle number in this case) in the wave 
function @,|0). 

For an even number of electrons the j/-pair wave func- 
tion in the projected form is 

l*'"^ = 2^ / "^^ ^"''"' n ( 1 + ^5k a|,X_ki ) |0) , 
•' k 

(21) 

where the contour is any counterclockwise path that en- 
closes the origin. If we perform the residue integral 
above, we recover the desired wave function Eq. (0). 

As written above the wave function is not normalized, 
so expectation values of observables will include a factor 
given by (^^|^^) in the denominator. A straightforward 
evaluation of this factor gives (^t,|^^ 
is defined by the residue integral |12[| 



i?[j, where RI^ 



i?™(ki,k2,---,k™) 



1 

27ri 



k^ki,k2,' 



\^ + i9l) (22) 



EE-- E 

Pl < P2 < <Pi/-7 



2 2 2 

5pi5p2---3p„ 



(23) 



where p^ ^ ki, k2, • • • , k^ (i = 1, • • • , i^ - n) 
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The integral in Eq. (g2|) has sharp oscillations as a func- 
tion of ^ (or 9 defined by the coordinate transformation 
^ = e'^) whose severity increases with increasing system 
size. We therefore evaluated the integral using the ana- 
lytical result Eq. (p3). This was accomplished efficiently 
with an algorithm that took advantage of the polynomial 
structure of the original integrand in (p2) , and led to a 
considerable speed up in the integral evaluations. 

These residue integrals are useful not only because the 
matrix elements can be easily evaluated, but also because 
they satisfy various recursion and sum formulas. For ex- 
ample, Dietrich et al. [O found the following recursion 
relation 



i?;7(ki,k2,---,k„)=i?™+^(ki,k2,---,k„,k) 

Tl-f 

1 + 



5^ii!:r++i'(ki,k2,---,k„,k). (24) 



This result follows straightforwardly from the definition. 
Another useful relation involves the derivative required 
in the variational principle: 



9Jt';r(ki,k2,---,k„ 
dgw 



= 2gki?:r+i(ki,k2,---,k,„,k) 



(25) 



In addition, we have also found and exploited the follow- 
ing sum rule: 

^ g2^:r+'i'(ki,k2,---,k„,k) 

k5^ki,k2,--,k„ 

= (i.-n)i?;r(ki,k2,---,k„) {i,>n). (26) 

All these relations follow in a very straightforward fash- 
ion from the definition. For completeness we include the 
case where the residue with two or more equal indices is 
required; it is then useful to rewrite the definition of the 
residue, Eq. (p2|) as 



(27) 



i?™(ki,k2,---,k™) 



2^^* J nk=ki,k.,.,k„ 



(1+^5^; 



The advantage of this definition over Eq. (^2|) is that 
there is no ambiguity when one or more momentum in- 
dex is used more than once. With this definition these 
residues still satisfy the recursion relation (|24|). This 
fact is useful for simplifying the expressions that involve 
restricted double momentum sums, which arise in the 
canonical formulation. 

Omitting details, the ground state energy can be writ- 
ten in a form reminiscent of the grand canonical BCS 
formulation (see Eq. (pi)): 



E, 



E2ek^ 

k 



9l 



■9i 



^i(k) 



N ^ 
k,k' 

-1 


ffk .9k' 


9l gl 



^?(k,k') 



r '2 



(k,k') 



■ 1 + 5^ 1 + 5^' 
where we have now defined normalized residues: 



(28) 



•^(l,,,...,],^)^^^^-2M Yl (1 + ^2,). 



i?g 



k'=ki,--,k„ 



(29) 



This definition is clearly motivated by the fact that in the 
bulk limit, the canonical results converge to the grand 
canonical ones. 

A straightforward but tedious variation of Eq. ( [2^ ) 
yields the following variational equation: 



( 2e\ -h Ak )5k = Ak [ 1 - 5^ 



where 



\U\ 



A 



N if i + gl '^ ' 



\u\ 



5k' 



(30) 

(31) 
(32) 



Ak^2ek5^}(k)[l-ri(k)] 



M 

N 



rl{k)il-gl)-rl{k,k) 



{l+gD Y, 2e^>-^Uik',k)-rlik)rlik')) 
k'#k ^ + 5k' V / 



\U\ 



k'^^^k 



9l 



2^5k E Yf^(r?(k',k)+gk'.gkr2^k',k) 



\U\ 



^(1 + 5^) E T 



gp gp' 

+ g^ 1 + .9 



r2'(p,p',k) 



'gpgp"'3(P:P'.k) 



M 

N 



{l+gl)rl{k)J2- 



gp gp' 
p,p. +5^1 + 5^' 

+.9pgp"'i(p,p') 



^l(p,p') 



(33) 



A solution to Eq. ( |30[) is obtained by numerical itera- 
tion. Note that Eq. (pQ) resembles the analogous equa- 
tion in the grand canonical ensemble, Eq. (ttSh. The fac- 
tor Ak appears here in addition; all terms in Eq. (p3) 
are of order 1/A^, and therefore vanish in the thermody- 
namic limit. Moreover, both the single-particle energy 
Eq. ( pl| ) and the pairing potential Eq. (B3) are modified 
by the normalized residue integrals. In practice, instead 
of Eqs. (p0|)-(|33|), we used the simpler- looking expression: 
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N denoin 

denom - (e^u - 2ek + ^"j i?}(k) 

iV 



E E 

P^kpVk, p'5^p 



.9p5p'^2(k,P,p') 



5^5p'^3(k,P,p') 



(34) 



Although this equation does not resemble the grand 
canonical equations, it is homogeneous in {^k} and has 
far fewer terms. 

In the case of an odd electron number N^ — 2iy+l, we 
define the fixed Ng wave function in terms of a residue 
integral as 



1^ 



21/+1/ 



1 

27ri 



d^C 



-'<. 



n 



(35) 



This wave function carries a momentum label q which 
gives the momentum of the unpaired electron (and hence 
of the total state). The normalization factor is now 
{^2iy+i\^2^+i) = -Ro(q), and the total energy is given 
by 



E2U+I — Eq 



\u\ 



E 

k/q 



2ek- 



\U\ 



9l 



r?(k,q) 






N '^ '^ 

k^^qkVq 



' TV 


y i + 5k ''o(q) 




5k 


.9k' 


r?(k,k', 


q) 


l + .9^ 


l + .9^' 


^^(q) 




9l 


9l 


r3(k,k', 


s) 



1+5^1 



9l 



rli^) 



(36) 



This expression for the ground state energy has the same 
form as Eq. (Pq ) except that the kinetic energy of the 
unpaired electron is singled out and the momentum sums 
explicitly prohibit the singly occupied momentum state. 
As in the even case a lengthy variational equation for the 
(7k 's is obtained and must be solved numerically. The 
end result is that we have the ground state energy for 
any number of electrons. 

Thus one can construct the gap by 



2A(iVe) = En^^i - 2En^^ + En^+1. 



(37) 



Various definitions of a gap or binding energy exist in 
the literature p^ , p^ ,p|,pl[ . The basic idea is the same 
— one wants to compare the difference in energies be- 
tween two systems, one in which 2iVe electrons are dis- 
tributed equally over two subsystems containing iVe elec- 
trons each, and the other in which the subsystems contain 



A^e + 1 and A^e — 1 electrons respectively. If A^e is even, 
the former has lower energy since pairing is fully utilized; 
this is reflected in a positive gap. If A^e is odd the lat- 
ter has lower energy and therefore the gap is negative. 
This is the origin of the so-called parity effect measured 
by Tinkham and coworkers |H-0,16l. For evaluating the 
gap A(A'e) in Eq. (p^), we choose the momentum q that 
yields the lowest energy for a given odd number of elec- 
trons. This is important for properly recovering the exact 
results in certain regimes, i.e. strong coupling and/or low 
electron density. 

In the grand canonical BCS formulation, there is a di- 
rect correspondence between the variational parameters 
and the occupation probability rik = X)cr('^k(T'^kcr)- It is 
determined through the gk's as 



5k 1 <=k 

rik = 2 - — !^ = I - 



9l 



E], 



(38) 



where the various functions are defined in Section II B 



Utilizing the Wk's and Vk's instead of the gk's makes the 
correspondence even more transparent, for we have rik = 
2f^ in that case. In the canonical formulation, for an 
even number of electrons, for example, we construct the 
matrix element 



?^kcr = (V'2iy| aLflktr \lp2iy) 



and obtain 



nk = 5]nk. = 2-f^ri(k) 



.9k 



(39) 



(40) 



2. Exact Solutions 

The formalism developed so far is applicable in any 
dimension. The results to be discussed later in this pa- 
per focus on one dimension only. One reason is that the 
model is very local, so that for the properties we will dis- 
cuss here, dimensionality is not too important. The sec- 
ond reason is that comparisons can be made with exact 
results, which are readily available only in one dimension. 

Ground state energies can be obtained both by ex- 
act diagonalization (on small system sizes) and by Bethe 
Ansatz techniques [[7|jlqjq|. We have already outlined in 
some detail B the numerical procedure used to obtain 
ground state energies for the attractive Hubbard model, 
and the reader is referred to those references for further 
details. 



III. RESULTS 
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A. Ground state energy 

We present results in one dimension and with only 
the nearest neighbours included for the electron hopping 
{ts =t). It has been shown in Ref. [g[ that for large sys- 
tems, the grand canonical BCS approximation yields the 
exact energy in the strong- and weak-coupling limits, and 
in the dilute limit for all coupling strengths. Thus devia- 
tions from the exact results are largest for weak to inter- 
mediate coupling strengths and for larger electron density 
n = {Ne)/N. We therefore anticipate improvements by 
the canonical BCS approximation for these cases, espe- 
cially for small system size N. Because of the particle- 
hole symmetry l7|q|, we need to study only up to half 
filling, n = 1 (i.e., the number of pairs is half the number 
of sites) . Note that for the exact and canonical calcula- 
tions, the density is defined simply by n = N^/N for a 
given number of electrons Ne- In the results shown below, 
for changing the density, we vary the electron number for 
a fixed system size, as is the case in the experiments. 

In Fig. 1, we show the ground state energy per site 
(absolute value) as a function of the electron density for 
(a) \U\/t = 10 and 4 and A^ = 16, and for (b) \U\/t = 2 
and A^ = 4, 8 and 32. The exact and canonical results 
are plotted with symbols and the grand canonical ones 
are shown with curves. In Fig. |l|(a), for the exact and 
canonical cases, the energies with even and odd num- 
bers of electrons are shown in different symbols. The 
conventional grand canonical BCS wave function con- 
tains only the components with even numbers of elec- 
trons. Strictly speaking, the exact and canonical results 
for odd numbers of electrons should be compared with 
the grand canonical ones in the parity-conserved scheme 
p8| , p3[ with the odd number parity. 

In Fig. 0(a), we first note the difference between the en- 
ergies with even and odd numbers of electrons, as clearly 
seen for \U\/t = 10. The difference becomes more ap- 
parent for larger coupling strengths, giving rise to the 
even-odd oscillations in the energy as a function of the 
electron density. For \U\/t = 10, we note another dif- 
ference between the even and odd electron numbers; the 
canonical energies are much closer to the exact ones for 
the even numbers. It can also be seen for this strong 
coupling case that the canonical results (for even A^e) 
are converged to the grand canonical ones for almost all 
density values, and both results are in very good agree- 
ment with the exact solutions for smaller density. As the 
coupling strength becomes smaller, for a fixed system 
size, the canonical energy deviates more from the grand 
canonical one and, for an even number of electrons, im- 
proves slightly the agreement with the exact energy for 
larger density. This can be seen for \U\/t = 4 in Fig. |](a), 
while the even-odd difference is now smaller. 
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FIG. 1. (a) Ground state energy as a function of the elec- 
tron density n = Ne/N for 16 sites. The points are the canon- 
ical and exact results (different symbols for even and odd A'^^), 
while the grand canonical ones are shown with curves. The 
difference between the energies with even and odd A^e can be 
seen clearly for \U\/t — 10. The canonical BCS results for 
even Ne axe better than those for odd Ne, and converge to 
the grand canonical energies for strong coupling. 
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FIG. 1. (b) Same as (a), but for j[/|/t = 2 and for N = 4,8 
(upper frame) and 32 (lower frame), and here the even and 
odd points are not distinguished. The improvement by the 
canonical scheme is more apparent for weak coupling and 
small system size, as can be seen for A = 4 and 8. For 
N — 32, the canonical and grand canonical results are more 
or less converged for all densities. 
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For smaller system size and coupling strength, we see 
more improvements due to the use of the canonical for- 
mulation. In the upper part of Fig. 0(b), we show the 
results for \U\/t = 2 and for iV = 4 and 8. Note that the 
energy range is magnified compared with Fig. 0(a), and 
that for the exact and canonical results, the even and 
odd energies are not distinguished by different symbols. 
For iV = 4 the agreement between the canonical and ex- 
act energies is excellent for all values of n, and it is still 
very good for N = 8. In the latter case, the grand canon- 
ical curve happens to be very close to the exact results for 
odd numbers of electrons for larger density. As explained 
above, however, the grand canonical results shown here 
must be compared for the even numbers only. In fact, the 
grand canonical energy for an even number of electrons 
differs most from the exact and canonical energies for 
larger density. On the other hand, the canonical energy 
converges to the grand canonical one as the system size 
becomes larger, as can be seen in Fig. |^(b) for N — 32. 
The even-odd difference is negligible for this weak cou- 
pling and large system size. 
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FIG. 2. (a) Ground state energy as a function of the cou- 
pling strength |f/|/f for N = 8 (upper frame) and 64 (lower 
frame) and for various values of n = N^/N . For N = 8, the 
improvement by the canonical method can be seen, and for 
n = 0.25, the canonical and exact energies are indistinguish- 
able. On the other hand, for A^ = 64, the canonical energies 
are converged to the grand canonical ones. Both the canonical 
and grand canonical results reproduce well the exact solutions 
for small n and in the zero- and strong-coupling limit. 

In Fig. 2, the ground state energy is plotted as a func- 
tion of the coupling strength for N = 8 and 64, for vari- 



ous densities with (a) even and (b) odd numbers of elec- 
trons. For small systems, the canonical results improve 
the grand canonical ones, especially for intermediate cou- 
pling strengths, while the former converge to the latter 
as the coupling strength is increased. This can be seen 
for TV = 8 in Fig. |(a). In fact for n = 0.25, the ex- 
act (solid) and canonical (dashed) lines are indistinguish- 
able, though the grand canonical one deviates from them 
only slightly. The three results converge as the coupling 
goes to zero, and also in the strong-coupling limit. For 
N = 64, the size is large enough that the canonical and 
grand canonical results are converged in the given scale 
for all densities and coupling strengths. 
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FIG. 2. (b) Same as (a), but here we compare the exact and 
canonical BCS results for an odd number of electrons A^e- For 
N — 8, the energies with odd A^e's are poorly reproduced for 
large \U\ by the canonical BCS approximation, compared with 
the even-A^e case shown in (a). As the system size is increased, 
the difference between the errors in the energy for even and 
odd A^e's becomes smaller, as can be seen for A'^ — 64. 

In Fig. H(b), the canonical and exact energies are com- 
pared for odd electron numbers. For N = 8, compared 
with Fig. 0(a), the exact energy is reproduced rather 
poorly by the canonical BCS approximation, especially 
for larger coupling strengths. Moreover, unlike the even 
number case, the canonical BCS results do not converge 
to the exact ones as the coupling strength is increased 
further. This difference between the errors in the energy 
by the even- and odd-iVg canonical BCS approximation 
turns out to be smaller as the system size becomes larger. 
For N — 64, the difference between the even and odd 
cases has diminished significantly from the TV = 8 case. 



B. Energy gap 
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1. Finite size effects 

We show the energy gap as a function of the density 
for (a) \U\/t = 1, (b) \U\/t = 1.5 and (c) \U\/t = 4; 
for A^ = 4 and 8 (upper figures) and iV = 16 (lower fig- 
ures) in Fig. H and for A^ = 32 (upper) and 64 (lower) in 
Fig. 0. Again, the exact and canonical results (A(A''e) in 
Eq. (|3^)) are plotted with symbols and the grand canoni- 
cal results (Ao in Eq. (16)) are shown with curves. In the 
upper parts of Fig. g, the exact solutions for A = 4 and 
8 are the circles and squares, respectively, and the corre- 
sponding canonical results are the crosses and stars. The 
gaps for odd numbers of electrons are negative. As dis- 
cussed above, the conventional BCS wave function that 
is used in this work does not contain the odd-A'e compo- 
nents, and thus the gap parameter defined by Eq. ( [iq ) 
misses out the gap for odd electron numbers. For even 
electron numbers, the canonical BCS method improves 
the grand canonical results significantly for weak cou- 
pling and small system size. 
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FIG. 3. (a) Energy gap as a function of the electron den- 
sity n = Ns/N for A^ — 4,8 (upper frame) and 16 (lower 
frame) and for \U\/t = 1. In the upper figure, the circles and 
crosses are the exact and canonical BCS results for A^ = 4, 
respectively; the squares and stars are for N = 8. The grand 
canonical results are shown with curves. The gap for odd 
Ne is negative, and for even Ne, the super- even oscillations 
{N^ — Am vs. 4m -I- 2) can be seen. The canonical results 
are in excellent agreement with the exact ones, while the 
grand canonical BCS result completely misses the gaps for 
Ne=4m + 2. 



FIG. 3. (b) Same as (a), but for \U\/t = 1.5. The symbols 
used for the exact and canonical gaps for A'' = 4 and 8 are 
the same as in (a). The canonical results are still in good 
agreement with the exact ones. 

In Fig. |(a) for |C/|/i = 1 and A^ = 4, 8 and 16, we 
can see the excellent agreement between the exact and 
canonical BCS results. Furthermore, there is a striking 
feature in these figures - the difference between the gaps 
for A^e — 4771 and Ag — Am + 2, where 777, is an integer. 
(Note that for all the results shown, the number of sites 
is a multiple of four.) The gaps for Ae — Am +2 are much 
larger than those for A'e — Am. Moreover, it is clearly 
seen for A^ = 16 that the gaps for A'e = 4777, -I- 2 increase 
as the density increases. We call these oscillations in the 
gap as a function of the even- A'e the super- even effect: 
the system is more stable with 4777 -f 2 electrons than 
with Am, electrons. This effect is more pronounced when 
the coupling is weak, and it stems from the quantized and 
doubly degenerate energy levels, Ck = e_k, of the unper- 
turbed system. As will be shown later, when the coupling 
is weak, the occupation probabilities of the unperturbed 
states can be approximated by those for zero coupling. 
We can thus understand the super-even effect as follows, 
in terms of the unperturbed energy levels which are oc- 
cupied by pairs of electrons up to the Fermi level and 
are empty above it. To simplify the discussion, we ignore 
the energy change due to the blocked states by unpaired 
electrons. 

In the case of one dimension and with ts = t, the 
unperturbed energy, Eq. (||), reduces to Cfc = —2t cos kR, 
where R is the lattice constant. From the periodic bound- 
ary condition, kR = 2n j/N {-N/2 < j < N/2), and 
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each energy level is degenerate for ±kR except for kR = 
and TT. Hence each level with < \kR\ < ir can accom- 
modate two pairs, (fc ti ~^ i) and {—k t, fc |). Thus 
when there are Am + 2 electrons, in the ground state, all 
the levels from kR = up to kpR are fully occupied by 
the pairs, while with 4m electrons, the Fermi level has 
a vacancy for one more pair. Therefore, when there are 
4m electrons, a way to break a pair with the minimum 
energy is to flip the spin of an electron of the pair at 
the Fermi level. The unpaired electrons then occupy the 
states ikp and there is no extra cost for the kinetic en- 
ergy. On the other hand, when there are 4m +2 electrons, 
one cannot break either pair at the Fermi level simply by 
flipping a spin, but an unpaired electron has to move up 
to the next available level, increasing the kinetic energy. 
This is why the gap for pair breaking is larger for Am + 2 
electrons than for 4m. electrons. The fact that the for- 
mer gap increases as a function of A^'e is particular to the 
one-dimensional band structure that we use: the level 
spacing becomes maximum around kR — 7r/2 (i.e., half 
filling) and so does the kinetic energy cost for breaking a 
pair. 
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FIG. 3. (c) Same as (a), but for \U\/t = 4. For N = 16 
the canonical gaps are converged to the grand canonical curve 
for almost all densities. Even for A'' = 4 and 8, the canon- 
ical results are closer to the grand canonical ones, and the 
super-even oscillations have disappeared. 

The super-even effect can be recognized more clearly in 
Fig. H(a) for A^ = 32 and 64. The overall scale of the gap 
becomes smaller for larger system size (note the reduced 
scale compared to Fig. 0(a)). For 32 sites, the canonical 



BCS results still follow the exact solutions closely for 
almost all density values. For 64 sites, the system is so 
large that even for this weak coupling, the canonical gaps 
are converged to the grand canonical curve in the dilute 
limit. The canonical results more or less follow the exact 
ones for n > 0.3, where the super-even effect is manifest. 



(a)|U|/t=1 
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FIG. 4. (a) Same as Fig. 3(a), but for larger system size; 
A'^ = 32 (upper frame) and A*' = 64 (lower frame). The su- 
per-even oscillations can be seen clearly, whereas the grand 
canonical BCS misses the 4m + 2 gaps completely. While for 
N=32 the canonical BCS results still follow the exact ones 
closely, for N=64 the former are closer to the grand canonical 
curve for low electron density. 

It can be seen in Figs. 3(a) and 4(a) that for larger 
density, where the 4m. and 4m. -I- 2 gaps are well sepa- 
rated, the grand canonical BCS completely misses the 
exact gaps for A^e = 4m -I- 2. On the other hand, the 
grand canonical curves appear to follow closely the ex- 
act 4m gaps. However, as the size is made larger than 
64 sites, the exact 4m gaps become smaller than the gap 
given by the grand canonical curve (see the bottom graph 
of Fig. ^a)) and finally in the bulk limit, they are quite 
small compared to the grand canonical values; at n — 1.0, 
the exact and grand canonical A/i are about 0.003 and 
0.015, respectively S. Meanwhile, the canonical 4m gaps 
go up towards the grand canonical curve, as they (and 
the 4m -I- 2 gaps) converge to the grand canonical values. 
Hence for some particular sizes (larger than 64 sites), 
the canonical results for the 4m gaps will be closer to the 
exact ones. 

In fact the grand canonical gaps tend to have a discon- 
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tinuity at iVg = Am + 2, but either value is generally well 
below the exact or canonical value. For smaller density 
where the super-even oscillations are not so prominent, 
the grand canonical gaps have cusps at N^; ~ Am -t- 2, as 
can be seen, e.g., for A^ = 64 and n < 0.5 in Fig. ^(a). 
Also, as the coupling strength is made a little larger but 
still small (|C/|/i ^ 2), these discontinuities at larger den- 
sity are replaced by cusps. This can be seen in Figs. 3(b) 
and 4(b) for \U\/t — 1.5. While for the small sizes shown 
in Fig. 0(b) there are still discontinuities for n > 0.5, for 
the larger sizes in Fig. 0(b) the curves are continuous for 
all densities, with cusps at N^ = Am, + 2. It is intriguing 
that the grand canonical BCS partially reproduces the 
super-even oscillations in this way. The grand canoni- 
cal solutions in the zero-coupling limit will be discussed 
further below. 
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FIG. 4. (b) Same as (a), but for \U\/t = 1.5. The canonical 
gaps are converged to the grand canonical ones for smaller 
density, while the former improve the latter still significantly 
near half filling. 

For A^ = 4 and 8 in Fig. 3(b), the agreement between 
the canonical and exact gaps is still very good, and the 
overall structure of the gaps as a function of the density 
is the same as for \U\/t = 1. For 16 sites, the canoni- 
cal results deviate slightly but still reproduce the exact 
gaps well. For small density, the canonical even-TVe gaps 
are closer to the grand canonical curve, while for larger 
density the canonical Ne = Am, + 2 gaps significantly im- 
prove the grand canonical curve. As the system becomes 
larger, the super-even oscillations diminish in amplitudes 
and also are confined more towards half filling. Also, as 



the size increases, for larger density where the (exact and 
canonical) Am, and Am + 2 gaps are well separated, the 
grand canonical curve shifts up relative to them, from 
near the Am gaps to above the Am, + 2 gaps. This can 
be seen in Figs. 3(b) and 4(b) in increasing order in size. 
For 32 sites, the last two cusps in the grand canonical 
curve happen to be closer to the exact values than the 
canonical ones. However, the canonical gaps capture the 
correct behaviour of the super-even oscillations. Finally 
for 64 sites, the canonical gaps are converged to the grand 
canonical one for low density, while they are much better 
for higher density, especially for N^ — Am. 

As the coupling strength is increased further (for a 
fixed size), the scale of the gap increases as a whole for 
both even and odd N^ (thus the even-odd difference be- 
comes larger), while the Am, vs. Am, + 2 difference de- 
creases. For example, it can be seen in Figs. 3(c) and 
4(c) that \U\/t = A is strong enough for most of the sizes 
shown for the system to reach the "bulk" limit, where 
the canonical and grand canonical methods hardly differ 
from one another. Even for A^ = 4 and 8, the super-even 
structure is gone, and the canonical gaps are closer to the 
grand canonical ones. 
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FIG. 4. (c) Same as (a), but for \U\/t = 4. For the large 
sizes shown, this is strong enough coupling so that the canoni- 
cal and grand canonical results are converged for all densities, 
and the super-even oscillations have disappeared. 

To see the effect of the coupling strength on the energy 
gap, we plot in Fig. Ha) the gap as a function of the 
coupling strength for N — 8 (upper) and 64 (lower), for 
quarter (n = 0.5) and half filling. Interestingly, for N — 8 
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the exact gaps for quarter and half filling are almost the 
same for all \U\ values, whereas for A^ = 64 they are 
somewhat different. On the contrary, in the BCS picture 
(both canonical and grand canonical), the difference be- 
tween quarter and half filling in the strong-coupling limit 
is about the same for small {N — 8) and large {N — 64) 
sizes and much larger than the exact result. As for the 
difference between the canonical and grand canonical re- 
sults, we first note that A'e = 4m for both cases shown in 
Fig. ra(a). For N = 8, the canonical and grand canonical 
results are equally good for weak coupling {\U\/t ^ 2), 
whereas the canonical gaps improve the grand canoni- 
cal ones slightly for stronger coupling. For large systems 
and for very weak coupling, the grand canonical results 
are better than the canonical ones. This can be seen for 
A?^ = 64 for |C7|/i < 1 for half filling, as we have seen in 
Fig. ||(a). As the coupling becomes stronger, the canon- 
ical gap converges to the grand canonical one, and this 
happens faster for larger systems. Indeed for 64 sites, the 
two curves for \U\/t > 2 can barely be distinguished in 
the given scale both for quarter and half filling, although 
half filling is a special case where the canonical gap de- 
fined by Eq. ( |37| ) does not converge to the conventional 
BCS gap - this will be discussed shortly. 





In Fig. ^(b), we compare the magnitude of the exact 
and canonical gaps for even and odd A^e for half filling, 
for 8 and 64 sites. For smaller systems, the magnitude of 
the gap for even N^ {— 4m) is a little larger than the one 
for odd Ne {— 4to ±1), as can be seen for N = 8. This 
difference is slightly larger for the exact solutions than 
the canonical BCS results. As the coupling strength or 
the system size is increased, this difference between the 
magnitude of the even (4to) and odd gaps diminishes. 
This can be seen in Fig. 0(b). Also note that the differ- 
ence between the BCS gap and the exact gap increases as 
the system size increases. (This can be seen in Fig. |^(a) 
as well.) This is due to the fact that the exact gap be- 
comes smaller for larger systems, whereas the canonical 
gap in the strong-coupling limit hardly changes as the 
system becomes larger (and A^ = 64 is large enough to 
be the bulk limit for \U\/t>A). 
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FIG. 5. (b) Same as (a), but here the magnitude of the 
gaps for even and odd electron numbers are compared for the 
exact and canonical BCS results. The comparison is made at 
half filling, so that even Ne = N and odd iVe = iV - 1. For 
N = 8 the even gaps are slightly larger in magnitude than the 
odd ones, while for A'^ = 64 the difference has diminished. 



FIG. 5. (a) Energy gap as a function of the coupling 
strength \U\/t for N = 8 (upper frame) and 64 (lower frame), 
for the densities n = 0.5 and 1.0. For N — 8 (and for 
Ne — 4m) the canonical and grand canonical results are 
equally good for weak coupling, while the former improve the 
latter slightly for stronger coupling. For N = 64 the canon- 
ical and grand canonical results can hardly be distinguished 
for almost all the coupling strengths shown. 



2. Grand canonical gap for weak coupling 

As mentioned above, for very weak coupling, the grand 
canonical gaps have discontinuities for A^e — Am -f 2 at 
larger density. We can understand how these discontinu- 
ities arise, by looking at the density as a function of the 
chemical potential in the zero-coupling limit. In Fig. ^(a) 
we show the density (top) and the gap Ao (middle) as a 
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function of the chemical potential fj,, and Ao as a function 
of the density n (bottom), for zero coupling and TV = 16. 
When \U\ = 0, Abcs -- 0, /l = ^i and Ek = |efe - A^l- 
There is no gap equation and n is simply determined by 
Eq. ( p^ for a given ^. Thus the existence of a gap is 
solely due to the finite system size. It can be seen in the 
top figure that as /i changes continuously, the density n 
changes as a step function: it is multi-valued when fj. is 
equal to any of the discrete Cfc due to the orbital and spin 
degeneracy. As n moves from one level to the next, the 
density stays the same, whereas Ao = nim{Ek) increases 
from zero, peaks when /z is precisely between the two lev- 
els, and falls to zero again. Hence Ao as a function of n 
has i5-function-like peaks, as seen in the bottom figure. 
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FIG. 6. (a) The electron density n (top frame) and the gap 
Ao (middle frame) as a function of the chemical potential /x, 
and Ao as a function of n (bottom frame), for A'^ = 16 and 
for zero coupling. The fact that the density is a step function 
of fj. and that the gap exists is because of quantized energy 
levels due to the finite system size. The Ao has peaks when 
|efc — /i| becomes a maximum. 



When the coupling is nonzero, n and fi must satisfy 
not only the number equation ( |l5|) but the gap equation 

1 1 -sr^ 1 

N ^ ~ 



- = -y 

\U\ /v Z^ 



2^k N 



-E 



2 v/(ek - lir + A^cs 



(41) 



which is obtained by substituting Eq. (nm into Eq. (12). 
It turns out that when the coupling is weak, fl given by 
Eq. (Ql]) cannot take certain values in between two levels 
(corresponding to the plateau regions in n in Fig. ra(a)), 
and this can be understood in the following way. 



For weak coupling, we can consider the relation be- 
tween the density and chemical potential (p.) roughly in 
the same way as for the zero-coupling case: for A'e = 4?7i, 
p. is very close to one of the levels efc (when \U\ = 0, it 
is equal to efc), while for N^ = 4to -|- 2 it must be in be- 
tween two levels. Therefore for N^ = Am, ek — p — and 
thus Abcs must be finite for n to have a certain value 
(see Eq. (p^). For iVg ~ Am + 2, however, if efc — /i 
is finite, then Abcs is driven to zero so that the RHS 
of Eq. (ETl) attains a large enough value for small \U\. 
Indeed, for sufficiently small |C/|, even if Abcs is zero, 
Eq. (Ell) cannot be satisfied for ft, near the middle of the 
range between two levels. Instead jl is driven towards 
either level, so that the RHS of Eq. (41) increases suffi- 
ciently to equal the LHS. The number equation (15) is 
still satisfied, since Abcs — 0. Thus certain values of jj. 
are not allowed for small values of \U\. This situation is 
illustrated in Fig. ^(b) for |{7|/i = 1. In the top part, the 
abrupt jumps in n seen in Fig. |^(a) have been somewhat 
smoothed, whereas the plateau parts have disappeared. 
Accordingly the gap becomes discontinuous as a function 
of jl as well as n. In the latter case, Ao indeed has two 
values for a given n at each discontinuity shown, corre- 
sponding to the two solutions for jl for coming up from a 
lower level and coming down from the next level. In any 
case, all values of n are possible (in contrast to jl). 



(b)N.16:|U|/t=1 




FIG. 6. (b) Same as (a), but for \U\/t — 1 and now the 
"chemical potential" is jl given by Eq. ([ll|). The jl cannot 
take certain values in between unperturbed energy levels, cor- 
responding to the plateau regions in (a). 

For smaller density or larger system size, the level spac- 
ings become smaller, and it becomes possible for jl to 
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have all the values in between levels. This also occurs 
for stronger coupling, as can be seen in Fig. @(c) for 
\U\/t = 2, where Aq is now continuous and has cusps 
at Nf. = 4m + 2. These cusps come about when the 
coupling is still weak enough so that Abcs is small, and 
for the same reason as in the zero-coupling case, that is, 
when \ek — Al becomes a maximum. 
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|k| carried by an electron and a hole are the same. In the 
canonical BCS scheme, for even TVe, the minimum pair- 
breaking energy is given by 2A(iVe) in Eq. (p7|), where 
the lowest energy is chosen for each of the systems with 
A^e — 1, -^e, and iVg -I- 1 electrons. For the lowest Sat^+i 
and i?Ar^_i, the momentum carried by an unpaired elec- 
tron and that by a hole, respectively, are not necessarily 
the same. (This is also the case for the exact solutions 
by the Bethe ansatz for a finite system.) In this section, 
we evaluate the gap A(iVe) by taking the same momen- 
tum for an electron and a hole as in the grand canonical 
picture, and compare it with the grand canonical quasi- 
particle energy. We should remark that the difference in 
gaps defined in these two different ways is generally quite 
small, but noticeable in certain extreme limits. 

In Fig. 1^, we plot A(7Ve) obtained by the canonical 
BCS method, for which the same \k\ has been taken for 
an electron and a hole and which we call A^ , as a function 
of the kinetic energy e^ . The energy band in one dimen- 
sion is from —It to 2t, corresponding to < |fci?| < tt. 
We show results for iV = 32 and for quarter and half fill- 
ing, i.e., iVe = 16 and Ae = 32, respectively; and for (a) 
\V\lt = 1, (b) \U\lt = 5 and (c) |C/|/t = 50. The canoni- 
cal results are shown with circles, while the grand canoni- 
cal energy E-^ = \/(Sk — A)^ + ^bcs ^^ plotted with solid 
curves. 



(a) N.32: lU|;t = 1 



FIG. 6. (c) Same as (a), but for |t/|/t = 2. This coupling 
strength is large enough so that p, can take all the values 
in between levels. On the other hand, the coupling is weak 
enough so that Abcs is small, and Ao has cusps when |efc — /I 
becomes a maximum, similarly to the zero-coupling case. 

It is intriguing that the weak but nonzero coupling pic- 
ture described above does not interpolate smoothly to the 
zero-coupling case. As |t/| decreases, the slope of n (cen- 
tred around A'e = 4m) as a function of \jl seen in Fig. 0(b) 
(top frame) will become sharper, while the correspond- 
ing curves in Ao iji) will be reduced to points. At exactly 
\U\ = 0, n becomes vertical around A'^e = \m for all m. 
However, all values of /i in between are allowed now, and 
Ao(/i) becomes continuous as in Fig. o(a) (middle frame). 



3. Quasi-particle energy 

In the conventional BCS theory, there is no distinc- 
tion between a particle- and a hole-excitation in that the 
energy E^ 



^Bcs ^^ ^^"^ same for either 



excitation |11|. As defined by Eq. (p^), the gap Ao is the 
lowest quasi-particle energy E^. The minimum energy re- 
quired for breaking a pair is 2Ao and thus, the momenta 




FIG. 7. (a) The momentum-dependent canonical gap Afc, 
for which the same |fc| is taken for an electron and a hole 
(circles), and the grand canonical quasi-particle energy S^ 
(curves), as a function of the kinetic energy et- The re- 
sults are for A'^ = 32 and \U\/t = 1, and for quarter (upper 
frame) and half (lower frame) filling. For this weak coupling, 
Ek ^ \ek — /i|. 
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We have seen in Fig. ^(a) that for \U\/t = 1 and n > 
0.5, the grand canonical gap (Ao) happens to be very 
close to the exact one for Nf. — Am, and for N — 32, 
also to the canonical gap (note that N^ = Am for both 
quarter and half filling). We also showed in Fig. H(c) 
that A^ = 32 is large enough for the canonical gap to 
converge to the grand canonical BCS gap for almost the 
entire density range, already for \U\/t = A. We find the 
result that the momentum-dependent canonical gap A^ 
defined in this section agrees remarkably well with the 
quasi-particle energy E^ (defined in the grand canonical 
context) for all quasi-particle m,om,enta (i.e., not just at 
the minimum in Figs. 7(a)-(c)). 

In the canonical picture, when the coupling is weak and 
there are Am, electrons, the lowest energy for breaking a 
pair is to create an electron and a hole both at the Fermi 
level. Thus |fc| for an electron and a hole is the same, 
and this explains the super-even effect. The gap with 
this configuration is the minimum value seen in Fig. |3(a) 
for both quarter and half filling. In the grand canonical 
scheme, Abcs becomes very small for this weak coupling 
and hence Ek — \ek — fl\, where /2 turns out to be almost 
at the Fermi level (slightly above) for quarter filling and 
it is at the Fermi level for half filling (/i = 0) . 



This happens sooner (i.e., for smaller \U\) for smaller 
number of electrons: not only the Fermi momentum for 
zero coupling is closer to zero than those for larger iVg, 
but also the optimal momentum starts shifting at weaker 
coupling. In the case of quarter filling for iV = 32, the op- 
timal momenta for N^ — 15 and 17 are both \kR\ — tt/A 
for zero couphng, while for \U\/t = 5, the ground state 
energies Ei^ and En have different momentum depen- 
dence and have a minimum at \kR\ = 37r/16 and 7r/4, 
respectively. Yet if we take the same momentum for both 
Ne = 15 and 17, A^ agrees very well with Ek for all the 
|fc|'s, as seen in Fig. ^(b). The A^ has its minimum when 
i?i5 -I- En is the smallest; this occurs at |fci?| = 37r/16 
for \U\/t = 5. For half filling, the optimal momenta for 
both Ne = 31 and 33 are still \kR\ = tt/2 {eu = 0) for 
\U\/t = 5, as in the zero-coupling case. 

The preceding discussion illustrates that to obtain the 
grand canonical limit, one must in general choose the 
quasi-particle momentum (identical for the particle and 
hole cases) to minimize the sum of the two odd-electron 
energies. For intermediate to strong coupling this will 
not be given by the momentum expected from the non- 
interacting limit (in Ref. H this procedure worked be- 
cause they adopted a particle-hole symmetric model - 
see below). The last case (next paragraph) shows this 
for a very strong coupling example. 



(b)N=32:|U|/t.5 






(c)N=32:|U|/t = 50 
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FIG. 7. (b) Same as (a), but for |)7|/f = 5. 



As the coupling strength increases, the momentum of 
an unpaired electron (or a hole) that yields the lowest en- 
ergy for an odd A"e system shifts from the Fermi momen- 
tum in the zero-coupling limit towards zero momentum - 
it eventually reaches zero momentum, that is, the bottom 
of the non-interacting band, in the strong-coupling limit. 



FIG. 7. (c) Same as (a), but for |t/|/f = 50. 

For extremely strong coupling such as \U\lt = 50, the 
optimal momenta for N^ — 15 and 17 are both zero. In 
such a case, the ground state energy increases almost 
linearly as a function of e/c from fci? = to tt. This can 
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be seen in the upper part of Fig. 0(d), where E15 and 
£^17 are the squares (with the left axis) and the crosses 
(with the right axis), respectively. Accordingly, A^ is 
minimum at kR = and increases linearly as a function 
of efc, as seen in Fig. |^(c). In the grand canonical case, 
p, (related to fj, by Eq. (pl|)) for such strong coupling is a 
large negative value, and Ek has a minimum well below 
the bottom of the band. 



(d)N=32:|U|/t = 50 




FIG. 7. (d) Ground state energy for odd Ne as a function 
of the quasi-particle kinetic energy efc, for A'' = 32. In the 
upper figure, the energies for A^e ~ 15 and 17 are shown 
with the squares (with the left axis) and the crosses (with the 
right axis), respectively; in the lower figure, the energies for 
Ne — 31 and 33 are plotted accordingly. 

Half filling is a special case due to the particle-hole 
symmetry. In this case p. is zero for any coupling strength 
and hence the quasi-particle energy is always minimum 
at |fci?| = 7r/2. By the canonical variation, the ground 
state energy for A^e = 31 (the largest odd number for 
n < 1) has its minimum at Stt/S. This can be seen in the 
lower part of Fig. 0(d) (the squares, with the left axis), 
although the differences in energy for different |fc|'s are 
rather small. The E33 (the crosses, with the right axis) is 
the mirror image of i^ai about \kR\ = 7r/2 as a function of 
the momentum (for any \U\), as follows from the particle- 
hole symmetry relation [Q. Interestingly, the optimal 
momentum for N^ = 31 (and for A'e = 33) is still close to 
7r/2 for this coupling strength. Indeed, one requires even 
stronger coupling to have E31 that behaves like Ei^ in 
Fig. rod), while _E33 will decrease in a symmetric fashion 
linearly from kR — to tt. 



As the coupling increases, the system of any finite size 
approaches the "bulk" limit. We have seen above for 
n — 0.5 that in the strong-coupling limit, the energies 
for both of the odd systems (i.e., Ng = 15 and 17) have 
their minimum values with the momentum A; = 0. This 
is true for any finite size and any density smaller than 
unity. Thus the canonical A(A'e) with the lowest i?jVe+i 
and Ej^^-i (N^ even) converges to Ao, which is, for such 
strong coupling, given by -\/(emin — p)^ + ^Ics- This is 
not the case for half filling, however: the lowest E31 and 
£■33 will have kR = and tt, respectively, while both 
should have kR = tt/2 for the canonical gap to converge 
to Ao. 

For Am + 2 electrons, the agreement of A^ with Ef^ is 
not as good for all quasi-particle momenta for weak cou- 
pling. In the strong-coupling limit, however, A^ eventu- 
ally converges to Ek for all the momenta as we have seen 
above for N,, — Am. Finally, as explained above for the 
super-even effect, the minimum gap for Ne — Am, + 2 for 
weak coupling has different |fc| 's for the Ae -I- 1 and A"e — 1 
systems. On the contrary, the grand canonical gap has 
always the same momentum for an electron and a hole. 
This explains why it does not reproduce the exact gap for 
A'e = Am + 2 for such weak coupling, while the canonical 
one does. 



C. Occupation probability 

To conclude the discussion of our results, we show in 
Fig. 0(a) the variational parameter {gk} and (b) the oc- 
cupation probability {uk} as a function of k and the 
coupling strength \U\, for A^ = 8 and 64 and for half 
filling. In both figures, the points at discrete values of k 
(and discrete values of |C/|) are simply connected by lines. 
The smallest \U\ shown in these figures is 0.25 f, and for 
A^ = 64, the increment in \U\ has been taken finer than 
for N — 8 (while the k increment is naturally smaller 
for larger system size) . The results shown have been ob- 
tained by the canonical BCS formalism. However, also 
in the grand canonical scheme, both {gk} and {nk} have 
the same overall shape as a function of k and \U\: the 
actual values of {gk} rnay be different but only slightly, 
so that they look the same as those in Fig. ^(a) in the 
given scale. 

We can see in Fig. ^(a) that in the strong-coupling 
limit, {gk} as a function of k is almost flat (in the given 
scale) for a given \U\, and thus all the unperturbed levels 
are almost equally mixed H. On the other hand, in the 
zero-coupling limit, {gk} becomes a cosine function of k. 
Since the grand canonical BCS yields the same behaviour 
of {gk}, this can be understood by the simple expression 
of Eq. (p^. As \U\ approaches zero, Abcs goes to zero, 
and Ek can be expanded up to the leading order in Abcs • 
Then Eq. (|l^) reduces to 
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5k = 



A 



BCS 



Ck - /*| 1 



A^ 



2 (ek - fi) 



rA2 



- (fk - A) 



5k 



and 



2|ek-Ml 



A 



BCS 



(42) 



Thus in the Umit of |[/| -> and Abcs —" 0, 
Abcs 



(fk > A) 



(43) 



5k ^ --: (ek < m) 



Hence gk ^ for ek > /i, while for ek < Mi 5k diverges 
in the zero-coupling limit, and provides an image of ek 
(note that for half filling, /i = 0). 



(a) variational parameters 



maximum and minimum values at \U\/t = 10 are ap- 
proximately the same for both cases, and Uk is about 0.7 
at fc = and about 0.3 at the edges of the band. As 
the coupling is made stronger, the occupation probabil- 
ity (and gk) becomes almost equal for all the states, and 
Uk — 0.5 for all k in the case of half filling. 

(b) occupation probability 







FIG. 8. (b) Occupation probability {rik} as a function of 
k and |(7|/t, for N — 8 (upper figure) and 64 (lower figure) 
and for half filling. At \U\/t — 10, Uk is roughly 0.7 and 
0.3 at fc = and vr, respectively (here the lattice constant 
R = 1). In the weak-coupling limit, the distribution for the 
non-interacting electron gas is recovered. 



FIG. 8. (a) Variational parameter {gk} as a function of k 
and the coupling strength \U\/t, for N = 8 (upper figure) and 
64 (lower figure) and for half filling. The increment in \U\/t 
has been taken to be finer for TV = 64, while the increment in 
k is naturally smaller. 

Furthermore, the occupation probability clearly shows 
how the distribution over the unperturbed states changes 
as a function of the coupling strength. It can be seen 
in Fig. 0(b) that the distribution function for the non- 
interacting case is recovered for weak coupling; Uk = 1 
for ek < fj- and for ek > fj-, and n/j = 0.5 at the doubly 
degenerate Fermi level. In Fig. Wh) , not only the scale is 
reduced compared to Fig. H(a), but also the relative (av- 
erage) height of Uk at \U\/t = 10 against the one in the 
zero-coupling limit is larger. Thus Uk may appear to have 
more variation as a function of k in the strong-coupling 
limit than gk does. However, the difference between the 



IV. SUMMARY AND DISCUSSIONS 

We have formulated BCS theory for a canonical en- 
semble, following earlier work by Dietrich et al. |12| and 
Falicov and Proctto [Q, and very recently by Braun 
and von Delft Q. We have also generalized the linear 
formulation introduced in Rcf. Ilfl] for any system size. 
However, this method has proven to be numerically too 
intensive for "large" systems, and provided at best only 
marginal improvement over the nonlinear canonical for- 
mulation. 

In this work we have adopted a very definite model, the 
attractive Hubbard model, for various reasons. First, we 
wanted to have an exact solution with which to monitor 
the improvement over the grand canonical scheme. These 
are available for the Hubbard model in one dimension by 
the Bethe Ansatz technique, and so we have used these as 
a benchmark throughout this work. Second, we wanted 
to study a system which, by choice of parameters, could 
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easily span the weak coupling to strong coupling regime, 
as well as the low density to high density limits. In this 
way we have observed the crossover from the bulk to 
quantum limit for a variety of regimes. The attractive 
Hubbard model is no doubt the "minimal" model that 
accomplishes this. Finally, we wanted to use a model 
which could readily be generalized to a realistic case, so, 
for example, one could use a parametrized tight-binding 
model to fit the band structure for Al to better describe 
ultrasmall metallic Al grains. 

In this work, however, we focused on the first two 
points listed above, with the goal of establishing generic 
trends. The existence of the parity effect, as already de- 
termined by a parity-projected grand canonical scheme 
||18| , P3| and more recently by the canonical formulation 
of a uniformly-spaced-level model |9[| , emerges quite nat- 
urally in this model. However, in addition to the even- 
odd effect, we have found a super- even effect, with os- 
cillations in the superconducting gap occurring between 
even electron numbers which arc multiples of 4 (Am) 
and non- multiples of 4 (Am + 2). The magnitude of 
the gap variation for even numbers of electrons is, in 
some cases, comparable to the even-odd variation, i.e. 
|A(4m) - A(4m + 1)| w |A(4m) - A(4m + 2)|. Thus, 
such oscillations should be observable in the same kind of 
tunneling experiments for seeing the even-odd effect. We 
note, however, a key ingredient for the super-even effect 
to occur is the double degeneracy of levels, which in our 
case comes about from the simple equality, et — c-k- In a 
more general case, the degeneracy structure will be more 
complicated, and therefore oscillations will exist but may 
not be as simply periodic as a function of electron num- 
ber as is the case here. The super-even effect is also 
a result of quantized energy levels due to finite system 
size, and the effect will be stronger for smaller systems. 
For very weak coupling, the grand canonical BCS fails 
to reproduce the super-even effect, while the canonical 
scheme does; indeed, it yields very good agreement with 
the exact solutions. 

Finally, we note that the grand canonical BCS quasi- 
particle dispersion relation was beautifully reproduced 
by the canonical results, simply by varying the odd num- 
ber ground state momentum, with the proviso that the 
electron and hole momenta were the same. This simple 
correspondence is quite surprising. 
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APPENDIX: LINEAR CANONICAL VARIATION 

We summarize the linearized formulation of the canon- 
ical variation. In addition to the wave function for 
A^e = 2i^ defined by Eq. dlq), we define the wave function 
for odd number of electrons N^ = 21^+1, following ||] as 

i*2.+i) = 4. ^ E •■■ E c^(ki,k2,...,k.) 



n «lT«-ka |0) , 



(Al) 



where spin of the extra electron a can be either up or 
down. The extra electron blocks the state q from being 
occupied by pairs. There is no variational parameter for 
the blocked state: for the variational calculations, we 
always choose q that gives the lowest energy. Note that 
the C's in the above equation are different from those for 
the same number of pairs in Eq. (Rq) due to the missing 

q 

We now derive the variational equation using |^i,) in 
Eq. (|l^) for an even number of electrons Ne — 2v. The 
formulae below are slightly modified for an odd number 
of electrons Np = 2i^+ f in terms of the odd wave function 



(Al). We mention the difference for N^ — 2v + 1 in the 
end. 

The normalization factor for Y^2v) defined in Eq. (|lq) 



(*2..|*2^) 



E--EE--E 

ki < < k„ pi < < p,^ 

C*(ki,...,k,)C(pi, •••,?,) 



|a_k„iak„T • ■ ■ a-kiiflki 



piT -pii 



4.T«-P.i 



|0) 



= E---E iG(ki,---,Mr. (A2) 

ki < < k„ 

Similarly, the expectation value of the kinetic energy op- 
erator m H = T + V is 

(*2.|r|*2.) = (*2.| E^kaLaka|*2.) 
kCT 

= E---E |C(ki,...,k,)p2(ek, + --- + ekj. (A3) 

ki < < k„ 

Calculation of the potential energy term is more tedious. 
After some lengthy operator algebra it can be written in 
a general form for v pairs as 
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{^2.\V\^2.) 



-^ (*2^| Yl 4T"-k+u"-k'+uak'T 1*2,.) 



kk'l 



ic^l 



ki < < k„ 

+ J2 C^*(p,ki,k2,---,k._i) 

P<ki 

+ Y. C*(ki,P,k2,---,k.-i) 

ki<p<k2 

+ •■• + E C*(ki,k2,---,k.-i,p) 

k„-i<p 

+ Y C'*(p,ki,k2,---,ki,_2,k^) 

P<ki 

+ Y C*{-ki,p,k2,---,K-2,K) 

ki<p<k2 

+ •■• + E '^*('^i'^2,---,k^_2,k^,p) 

k„<p 

+ ••• 

+ J2 C*(p,k2,k3,---,k.) 

P<k2 

+ J2 C^*(k2,P,k3,---,k,) 

k2<p<k3 

+ •■• + E ^*(''2,k3,---,k„p)] (A4) 

k„<p 

The first term on the right hand side corresponds to a 
Hartree-like term, and gives the interaction energy of 
each of the i^ pairs with the other i/ — 1 pairs. The other 
terms systematically consider the various cases when one 
pair is scattered into another (unoccupied) pair state. 

In the following, we take the C's to be real variables. 
In Ref. [|2[ it has been proved for the BCS formulation 
that for negative pairing-type interactions, real (and pos- 
itive) variational parameters yield the lowest energy. The 
minimization condition for the energy ( po[ ) with respect 
to the real {C} results in the equation 

C(ki, • . . , k,) = -i M [ ^ C(p, ki, k2, • • • , k,_i) 

P<ki 

+ E C(ki,P,k2,---,k^_i) 

ki<p<k2 

+ ••• + E C(ki,k2,---,k,_i,p) 

k,.-i<p#k„ 

+ Y (^(p, ki,k2,---,k^_2,k^) 

P<ki 

+ Y C(ki,p,k2,---,k^_2,ki.) 

ki<p<k2 
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Y C(ki,k2,---,k^_2,P,k^) 

k„_2<P#k„_i<k„ 

Y C(ki,k2,---,ki._2,k^,p) 

k„<p 



E C'(p,k2,k3,---,k,,) 

p#ki<k2 

E C'(k2,p,k3,---,kj.) 

k2<p<k3 

Y C(k2,k3,---,ki.,p) 

k„<p 



where 



^. -2(ek, 



^k. 



M 

N 



(A5) 



(12') 



This should be solved for all the C's selfconsistently. 

For an odd number of electrons, in the energy ex- 
pectation value and the variational equation in terms of 
I'^iiy+i), the blocked q is excluded in all the {k} sums 
and {p} sums. Other than that, the normalization factor 
for |^2i/-i-i) is the same as Eg. (|A2|): the kinetic energy 
term is the same as Eq. ( |A3| ) except that 2(eki +• ' • + £k„) 
should be replaced by 2(eki + ■ • • + eO_+ Cq: the poten- 



tial energy term is the same as Eq. (A4) except that the 
factor ly {ly — 1) should be replaced by v^: the variational 
equation is the same as Eq. (A5), while in d as defined in 
Eq. ( |A6[ ) — Eq should be added and ly^ should be replaced 

We have found that this linear variation in terms of 
the C's just barely improves the ground state energy, 
compared to the nonlinear variation in terms of the g's. 
Interestingly, the C-formulation does not improve the 
^-formulation for the energy gap. We illustrate this in 
Fig. I for iV = 16 and for \U\/t = 1.5 and 4. We found 
that for smaller system size, the two formulations do not 
make any difference for any coupling strength, for the 
gap as well as the ground state energy. For larger sys- 
tem size (which was limited for N ^ 30) the gap from 
the C-formulation is slightly worse than the one from the 
f;-formulation for weak coupling and larger density. This 
can be seen in Fig. |9|(a), while for \U\/t = 4 in (b) the 
two results have converged for all the density. 
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X CBCS 4 LCBCS 
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(b) |U|/t = 4 

o o o o o o 

exact X CBCS » LCBCS 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



FIG. 9. Energy gap as a function of the electron density n 
for iV = 16, for (a) \U\/t = 1.5 and (b) \U\/t = 4. The exact 
solutions are shown with circles, while the canonical results 
from the g-formulation and from the linear C-formulation are 
shown with the crosses and triangles, respectively. 
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